#!/usr/bin/python
# GM_EstimateDependence.py
# Yet another module it looks like I never finished.  I was maybe trying to estimate the slope of that one big honking plot from Han, because ours didn't match?  I'm not sure.

import math
import numpy

def ED(simpars, phypars, appearance, vtuple):
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    (vrelappl, sigmaappl2, sigmaappb2) = vtuple
    # Distances
    # 1d
    cosb = numpy.cos(ba)
    cosb_db = cosb*deltaba
    d_d_dd = d**2*deltad
    d_d_d_dd = d**3*deltad
    # 2d
    dOmega = numpy.multiply(deltala[:,numpy.newaxis], cosb_db[numpy.newaxis,:])
    # 3d
    volelement = numpy.multiply(d_d_dd[:,numpy.newaxis,numpy.newaxis], dOmega[numpy.newaxis,:,:])
    dist_volelement = numpy.multiply(d_d_d_dd[:,numpy.newaxis,numpy.newaxis], dOmega[numpy.newaxis,:,:])
    # 0d
    Omega = dOmega.sum()
    del cosb, cosb_db, d_d_dd, d_d_d_dd
    # Integrate over locations of apparent magnitudes to obtain typical distance
    numberlist = []
    distancelist = []
    vtapplist = []
    proptolist = []
    for i in range(0, len(appearance)):
        # Find overall numbers for each apparent magnitude bin and typical distances thereof
        dist_number_integrand = dist_volelement*appearance[i]
        dist_number = dist_number_integrand.sum()
        number_integrand = volelement*appearance[i]
        number = number_integrand.sum()
        numberlist.append(number)
        dist = dist_number/number
        distancelist.append(dist)
        halfdist = 0.5*dist
        dindex = dist/PtScaled
        hdindex = halfdist/PtScaled
        # Use that knowledge to find typical relative angular speeds
        vtSintegrand = vrelappl[dindex,:,:]*dOmega
        vtS = vtSintegrand.sum()/Omega
        sig2lSgrand = sigmaappl2[dindex,:,:]*dOmega
        sig2lS = sig2lSgrand.sum()/Omega
        sig2bSgrand = sigmaappb2[dindex,:,:]*dOmega
        sig2bS = sig2bSgrand.sum()/Omega
        vtLintegrand = vrelappl[hdindex,:,:]*dOmega
        vtL = vtLintegrand.sum()/Omega
        sig2lLgrand = sigmaappl2[hdindex,:,:]*dOmega
        sig2lL = sig2lLgrand.sum()/Omega
        sig2bLgrand = sigmaappb2[hdindex,:,:]*dOmega
        sig2bL = sig2bLgrand.sum()/Omega
        vtapp = ((vtS - vtL)**2 + sig2lS + sig2bS + sig2lL + sig2bL)**0.5
        vtapplist.append(vtapp)
        propto = number*vtapp*dist**0.5
        proptolist.append(propto)
    numberlist = numpy.array(numberlist)
    distancelist = numpy.array(distancelist)
    vtapplist = numpy.array(vtapplist)
    proptolist = numpy.array(proptolist)
    return (numberlist, distancelist, vtapplist, proptolist)
